clc,clear;
close all;
warning off;
%导入数据
[~,name]=xlsread('newdata.xlsx');
data0=readmatrix('newdata.xlsx');
data0(isnan(data0))=10000;%为方便看
%%  导入数据(时间序列数据一列)
% for i=3:size(data0,2)
t=2;
data=data0(:,2+t);
%% 样本
nums = length(data);   %样本数
history = 5;    %历史特征个数
% number=input('输出数为：');%要与history成倍数
number=1;
%% 划分数据集
xd=zeros(1,history+number);
for i=1:nums-history-number+1
    xd=data(i:i+history-1)';
    for j=1:number
        xd=[xd,data(i+history+j-1)];
    end
    re(i,:)=xd;
end

%%  数据归一化 索引
X=re(:,1:end-number);
Y=re(:,end-number+1:end);
[x,psin]=mapminmax(X',0,1);   %归一化
[y,pcout]=mapminmax(Y',0,1);   %归一化

%% 划分训练集与测试集
num=length(re); %总样本数
% k=input('是否打乱样本(是1，不是0):');
k=1;
if k==0
    state=1:num;
else
    state=randperm(num);
end
r=0.7; %训练集占比
trainnum=floor(num*r);%训练集总样本数
xtrain=x(:,state(1:trainnum));
ytrain=y(:,state(1:trainnum))';

xtext=x(:,state(trainnum+1:end));
ytext=y(:,state(trainnum+1:end))';

%% 适应格式
for i=1:trainnum
    x_train{i,1}=xtrain(:,i);
end

for i=1:num-trainnum
    x_text{i,1}=xtext(:,i);
end

%% 创建模型
layers=[
    sequenceInputLayer(history)             %输入层
    lstmLayer(6,'OutputMode','last')        %LSTM层
    reluLayer                               %Relu激活层
    fullyConnectedLayer(number)             %全连接层
    regressionLayer ];                      %回归层

%% 参数设置
options=trainingOptions('adam' ,...                              %梯度下降算法
                        'MaxEpochs',500, ...                     %最大迭代次数
                        'InitialLearnRate',0.01, ...             %初始学习率
                        'LearnRateSchedule','piecewise', ...     %学习率下降
                        'LearnRateDropFactor',0.1,...            %学习率下降因子
                        'LearnRateDropPeriod',200, ...           %经过多少次训练后学习率为*0.1
                        'Shuffle','every-epoch', ...             %每次训练打乱数据集
                        'Plots','training-progress', ...         %画出曲线
                        'Verbose',true);                         %观察进度

%% 训练模型
net=trainNetwork(x_train,ytrain,layers,options);

%% 仿真预测
re1=predict(net,x_train);
re2=predict(net,x_text);

%% 反归一化（实际值）
Ytrain=Y(state(1:trainnum),:)';
Ytext=Y(state(trainnum+1:end),:)';

%预测值
pre1=mapminmax('reverse',re1',pcout);
pre2=mapminmax('reverse',re2',pcout);

%% 预测
futurenum=5; %预测数要和number成倍数
u=length(y)-history/number+1;
%按number
oxdi=y(1,u);
j=0;
for i=1:number
    if i<2
        continue;
    end
    oxdi=[oxdi,y(number,u)];
    j=j+1;
end
for i=u+1:length(y)
   oxdi=[oxdi,y(1,i)]; 
   for k=1:j
       if j==0
           break;
       end
       oxdi=[oxdi,y(k,i)]; 
   end
end
inputnew{1,:}=oxdi;%输入历史数据(进行下一步预测的)
outputs=zeros(futurenum/number,number); %初始化未来预测值
for i=1:futurenum/number
    %以历史数据作为输入
    indata=inputnew{i};
    outputs(i,1:number)=predict(net,indata(:));
    %更新历史数据
    inputnew{i+1,:}=[indata(1+number:end),outputs(i,1:end)];
end
%反归一化
reslout=mapminmax('reverse',outputs',pcout);

%% 未来预测图
u=futurenum/number;
oxdit=reslout(1:end,1)';
for i=2:u
   oxdit=[oxdit,reslout(1:end,i)']; 
end
figure
plot(1:nums,data,'r-o','LineWidth',1);
hold on
plot(nums:nums+futurenum,[data(end);oxdit'],'b-*','LineWidth',1);
set(gca,'XTick',1:nums+futurenum,'XTickLabel',[data0(:,1);(2024:2028)']);
xlabel('年份');
ylabel('医疗卫生机构数');
legend('历史数据','未来预测数据');
title([cell2mat(name(1,2+t)),'预测结果'])
grid on

%% 分别计算y
for i=1:number
        disp(['第',num2str(i),'个输出指标'])
        p1=pre1(i,:);
        p2=pre2(i,:);
        Ytr=Ytrain(i,:);
        Yte=Ytext(i,:);
        %% 均方根误差
        error1= sqrt(sum((p1-Ytr).^2)./trainnum);
        error2= sqrt(sum((p2-Yte).^2)./(num-trainnum));

        %% 相关指标计算
        %R^2
        R1=1-norm(Ytr-p1).^2/norm(Ytr-mean(Ytr)).^2;
        R2=1-norm(Yte-p2).^2/norm(Yte-mean(Yte)).^2;
        
        % MAE
        mae1=mean(abs(Ytr-p1));
        mae2=mean(abs(Yte-p2));
        
        %% 图
        figure 
        plot(1:trainnum,Ytr,'r-^',1:trainnum,p1,'b-^');
        legend('真实值','预测值');
        ylabel('医疗卫生机构数');
        title('训练集预测结果对比');
        
        figure 
        plot(1:num-trainnum,Yte,'r-^',1:num-trainnum,p2,'b-^');
        legend('真实值','预测值');
        ylabel('医疗卫生机构数');
        title('测试集预测结果对比');
        
        %百分比误差图
        figure 
        plot((p1-Ytr)./Ytr,'b-o','LineWidth',1);
        legend('百分比误差');
        ylabel('误差');
        title('训练集百分比误差曲线');
        
        figure 
        plot((p2-Yte)./Yte,'b-o','LineWidth',1);
        legend('百分比误差');
        ylabel('误差');
        title('测试集百分比误差曲线');
        
        %拟合图(相关系数r)
        figure
        plotregression(Ytr,p1,'训练集');
        set(gcf,'Toolbar','figure');
        hold on
        figure
        plotregression(Yte,p2,'测试集');
        set(gcf,'Toolbar','figure');
end
% end